###### Libraries ######
library(foreign)
library(arm)
library(Hmisc)
library(ggplot2)
library(memisc)
library(apsrtable)
library(car)
library(Amelia)

##### Data #####
#d <- read.csv('Desktop/Pesquisas/Legislative Behavior (FP1)/Data set/Lula03-10.csv', sep=',', header=T)
#d <- d[d$year==2005, ]
#d2 <- read.csv('Desktop/Pesquisas/Artigo Corrupção/bancocorrupcao.csv', sep=',', header=T)
#a <- merge(d2, d, by="name", all.x=T)
#write.csv(a, 'Desktop/Pesquisas/Artigo Corrupção/Revisao/bcorrup.csv')
#rm(a, d, d2)

#### Data ####
dat <- read.csv("corrupcao.csv", header=T, sep=",", na.strings="NA")
names(dat)

dat$escolaridade <- as.numeric(dat$GraudeInstruCAo)
dat$esco[dat$escolaridade==5] <- 0 
dat$esco[dat$escolaridade==2] <- 1
dat$esco[dat$escolaridade==1] <- 2
dat$esco[dat$escolaridade==4] <- 3
dat$esco[dat$escolaridade==3] <- 4
dat$esco[dat$escolaridade==7] <- 5
dat$esco[dat$escolaridade==6] <- 6

#### Efeito da Variação Estadual ####
fit1 <- glm(formula = reeleito ~ factor(UF), family=binomial(link=logit), data=dat) 
summary(fit1)

# Graph 1 #
est <- data.frame(estados=unique(dat$UF), regiao=c(0,1,0,0,1,1,4,2,4,1,2,4,4,0,1,1,1,3,2,1,0,0,3,3,1,2,0), coef(fit1), se.coef(fit1))
colnames(est)[3] <- "fit"; colnames(est)[4] <- "se.fit"
est <- transform(est, ylo= fit-(1.96*se.fit), yhi= fit+(1.96*se.fit))
est <- est[order(est$fit),]

pdf(file="estados.pdf", width=12, height=7)
ggplot(est, aes(x=as.factor(reorder(estados, rank(fit))), y=fit, ymin=ylo, ymax=yhi)) + geom_pointrange(size=1.5) + coord_flip() + xlab('Estados') + ylab('Efeito sobre Reeleição')+ labs(title="") + theme_bw() + theme(axis.title.x=element_text(face="bold", size = 18), axis.text.x=element_text(vjust=0.5, size=19)) + theme(axis.title.y=element_text(face="bold", size=16), axis.text.y=element_text(vjust=0.5, size=16)) + theme(legend.text=element_text(colour="black", size=13), legend.position="right")
dev.off()

coef(fit1)

#### Efeito Mensalao ####
fit2 <- glm(formula = reeleito ~ pobreza, family=binomial(link=probit), data=dat) 
summary(fit2)

fit3 <- glm(formula = reeleito ~ Sexo + Idade + esco + pobreza, family=binomial(link=probit), data=dat) 
summary(fit3)

fit4 <- glm(formula = reeleito ~ ideologia + aumevot0206 + Migrante + pork + Cadeiras_por_estado + despesa, family=binomial(link=probit), data=dat) 
summary(fit4)

fit5 <- glm(formula = reeleito ~ Sexo + Idade + esco + pobreza + ideologia + aumevot0206 + Migrante + pork + Cadeiras_por_estado + despesa + escandalo2, family=binomial(link=probit), data=dat) 
summary(fit5)

dat <- within(dat, Escandalo <- relevel(Escandalo, ref = 3))

fit6 <- glm(formula = reeleito ~ Sexo + Idade + esco + pobreza + ideologia + aumevot0206 + Migrante + pork + Cadeiras_por_estado + despesa + Escandalo, family=binomial(link=probit), data=dat) 
summary(fit6)

mtable(fit2,fit3,fit4,fit5,fit6)


##### Graphs ####
# Graph 2 #
a <- table(dat$Escandalo,dat$reeleito,dat$ideologia)

votos <- factor(rep(c('Esquerda','Centro','Direita'), c(3,3,3)),levels=c('Esquerda','Centro','Direita'))
escand <- ordered(rep(c('Não','Mensalão','Sanguessuga'),3),levels=c('Não','Mensalão','Sanguessuga'))
noreelt <- c(48,3,3,39,0,4,45,1,35)
reelt <- c(77,3,1,98,0,1,88,2,2)

logit.reelt <- log(reelt/noreelt)
logit.reelt[5] <- 0
data.frame(votos,escand,noreelt,reelt,logit=logit.reelt)

pdf(file="logit.pdf", width=12, height=7)
plot(rep(1:3, 3), c(-3.5,-3,-2.5,-2,-1.5,-1,-.5,0,1), type='n', axes=F, main='',xlab='Escândalo', ylab='Logit (Reel / Não-Reel)')
axis(1, at=1:3, labels=c('Nenhum','Mensalão','Sanguessuga')) # x-axis
axis(2, at=c(-4:4)) # y-axis
box()

points(1:3, logit.reelt[1:3], pch=1, type='b', lty=1, lwd=3, cex=2, col='black')
points(1:3, logit.reelt[4:6], pch=16, type='b', lty=2, lwd=3, cex=2, col='black')
points(1:3, logit.reelt[7:9], pch=8, type='b', lty=3, lwd=3, cex=2)
legend(c(-1), c("Esquerda","Centro","Direita"), lty=c(1,2,3), pch=c(1,16,8),lwd=c(3,3,3),col=c('black','black','black'))
dev.off()

# Graph 3 #
mod7 <- update(fit6, .~. - (Sexo + Idade + esco + Migrante + pork + despesa + Cadeiras_por_estado))

predictors <- expand.grid(pobreza=seq(.1,.7,.1),ideologia=0:2,aumevot0206=0:1,Escandalo=c("MensalAo","MensalAo e Sanguesuga","NAo","Sanguesuga"))
p.reel <- predict(mod7, predictors, type='response')
p.hat <- cbind(predictors, p.reel)

pdf(file="estima.pdf", width=6, height=10)
par(mfrow=c(4,1)) # 2 row and 2 columns
# Esquerda - não aumentou #
plot(c(.1,.7), c(0,1), type='n', xlab='Pobreza', ylab='Probabilidade Estimada', main='Esquerda - Não Aumentou Votação')
#lines(seq(.1,.7,.1), p.reel[85:91], lty=2, lwd=3, pch=4, type='b') # Nenhum
lines(seq(.1,.7,.1), p.reel[1:7], lty=1, lwd=3, pch=1, type='b') # Mensalao
lines(seq(.1,.7,.1), p.reel[127:133], lty=3, lwd=3, pch=16, type='b') # Sanguessuga
#lines(seq(.1,.7,.1), p.reel[43:49], lty=5, lwd=3, pch=3, type='b') # Mens+Sangues
legend("topleft", c("Mensalão","Sanguessugas"), lty=c(1,3), pch=c(1,16),lwd=c(3,3))

# Esquerda - aumentou #
plot(c(.1,.7), c(0,1), type='n', xlab='Pobreza', ylab='Probabilidade Estimada', main='Esquerda - Aumentou Votação')
#lines(seq(.1,.7,.1), p.reel[106:112], lty=2, lwd=3, pch=4, type='b') # Nenhum
lines(seq(.1,.7,.1), p.reel[22:28], lty=1, lwd=3, pch=1, type='b') # Mensalao
lines(seq(.1,.7,.1), p.reel[148:154], lty=3, lwd=3, pch=16, type='b') # Sanguessuga
#lines(seq(.1,.7,.1), p.reel[64:70], lty=5, lwd=3, pch=3, type='b') # Mens+Sangues

# Direita - não aumentou #
plot(c(.1,.7), c(0,1), type='n', xlab='Pobreza', ylab='Probabilidade Estimada', main='Direita - Não Aumentou Votação')
#lines(seq(.1,.7,.1), p.reel[99:105], lty=2, lwd=3, pch=4, type='b') # Nenhum
lines(seq(.1,.7,.1), p.reel[15:21], lty=1, lwd=3, pch=1, type='b') # Mensalao
lines(seq(.1,.7,.1), p.reel[141:147], lty=3, lwd=3, pch=16, type='b') # Sanguessuga
#lines(seq(.1,.7,.1), p.reel[57:63], lty=5, lwd=3, pch=3, type='b') # Mens+Sangues

# Direita - aumentou #
plot(c(.1,.7), c(0,1), type='n', xlab='Pobreza', ylab='Probabilidade Estimada', main='Direita - Aumentou Votação')
#lines(seq(.1,.7,.1), p.reel[120:126], lty=2, lwd=3, pch=4, type='b') # Nenhum
lines(seq(.1,.7,.1), p.reel[36:42], lty=1, lwd=3, pch=1, type='b') # Mensalao
lines(seq(.1,.7,.1), p.reel[162:168], lty=3, lwd=3, pch=16, type='b') # Sanguessuga
#lines(seq(.1,.7,.1), p.reel[78:84], lty=5, lwd=3, pch=3, type='b') # Mens+Sangues
dev.off()
